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Abstract 

The majority of nuclei available for study in solid state Nuclear Magnetic Resonance 
have half- integer spin I > 1/2, with corresponding electric quadrupole moment. As 
such, they may couple with a surrounding electric field gradient. This effect intro- 
duces anisotropic line broadening to spectra, arising from distinct chemical species 
within polycrystalline solids. In Multiple Quantum Magic Angle Spinning (MQ- 
MAS) experiments, a second frequency dimension is created, devoid of quadrupolar 
anisotropy. As a result, the center of gravity of peaks in the high resolution dimen- 
sion is a function of isotropic second order quadrupole and chemical shift alone. 
However, for complex materials, these parameters take on a stochastic nature due 
in turn to structural and chemical disorder. Lineshapes may still overlap in the 
isotropic dimension, complicating the task of assignment and interpretation. A dis- 
tributed computational approach is presented here which permits simulation of the 
two-dimensional MQMAS spectrum, generated by random variates from model dis- 
tributions of isotropic chemical and quadrupole shifts. Owing to the non-convex 
nature of the residual sum of squares (RSS) function between experimental and 
simulated spectra, simulated annealing is used to optimize the simulation param- 
eters. In this manner, local chemical environments for disordered materials may 
be characterized, and via a re-sampling approach, error estimates for parameters 
produced. 
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1 Introduction 



Since the discovery of Nuclear Magnetic Resonance (NMR), there has been 
great interest in the study of quadrupolar nuclei. These nuclei have an electric 
quadrupole moment Q, which couples with a non-zero electric field gradient. 
As a result, anisotropic frequency dependence is introduced, promoting over- 
lap between lineshapes arising from distinct chemical sites in powdered solids 
and degrading resolution. This issue has been addressed over the course of time 
by a number of experimental approaches. Early in the development of solid 
state NMR, Magic Angle Spinning (MAS) [1] was proposed, which reduces or 
eliminates second rank interaction terms and therefore broadening associated 
with the first order quadrupole interaction. This interaction depends explicitly 
on the angle 9 between sample rotor axis and the static, applied field of NMR. 
Attention here is restricted to first and second order quadrupole effects, each of 
which is a function of the second order Legendre polynomial Pzi®)- The second 
order quadrupole perturbation is also a function of the fourth order Legendre 
polynomial P±{6). Additionally, an appreciable second order isotropic shift in- 
versely proportional to the Larmor frequency ujq occurs; the center of gravity 
of a quadrupole lineshape is subsequently changed from the chemically shifted 
value. The characteristic features of quadrupole spectra provide valuable local 
bonding information and hence extensive work has been devoted to both re- 
solving individual chemical sites, as well as lineshape simulation. However, if 
the magnitude of the quadrupole interaction is significant, spinning sideband 
manifolds arising from satellite frequency transitions may still obscure spectra 
in one dimension [2]. Double Rotation (DOR) [3 J and Dynamic Angle Spinning 
(DAS) [Ifo] are successful in eliminating the effects of both second and fourth 
rank tensor terms, and thus also second order quadrupole broadening. More 
recently, Multiple Quantum Magic Angle Spinning (MQMAS) [6117] and Satel- 
lite Transition Magic Angle Spinning (STMAS) [81I91TT0] have become popular 
owing to mechanical simplicity. These procedures involve collecting data as a 
function of two independent time intervals in the pulse sequence [11] under 
Magic Angle Spinning conditions. Within the MQMAS experiment, directly 
observable single quantum coherence frequency transitions are correlated with 
multiple quantum transitions [12] and, in the case of STMAS, satellite transi- 
tions, which evolve between pulses and are selected via an appropriate phase 
cycle. From the center of gravity of peaks along the high resolution axis, 
isotropic shifts are deduced which are a function of both isotropic chemical 
(S l c s s °) and second order quadrupole (51£q) shifts. In turn, the isotropic second 
order quadrupole shift is a function of both the quadrupole coupling constant 
C q and asymmetry parameter r] q . The importance of these quantities lies in 
the fact that they are functions of the electric field gradient tensor V, and thus 
the details of the local bonding environment: 
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In order to unequivocally determine both C q and rj q , simulation of experi- 
mental spectra is necessary [13]. In the case of disordered chemical environ- 
ments [T^1I15II16|[TT] . calculations of powdered lineshapes for MQMAS becomes 
a formidable task. This is due to the fact that parameters relevant to simula- 
tion take on a distributed nature [TSITT^] . The focus of this paper is devoted to 
the optimized simulation of multiple quantum magic angle spinning spectra, 
in the presence of low to significant disorder. This is accomplished using quasi- 
random numbers sampled from model distributions of isotropic chemical shift 
and quadrupole coupling constant. Simulated annealing is used to optimize the 
non-convex RSS function, and in distinction to existing simulation methods, 
model parameter error estimates are calculated, using the non-linear jack- 
knife [20]. The overall process has been implemented in the C programming 
language with some tasks performed using the OCTAVE scripting language, 
and is highly amenable to distributed computing |21j. 



2 Theoretical background 

2.1 Lineshape Simulation 

Since the introduction of MQMAS experiments, there have been significant 
improvements in excitation efficiency and coherence transfer, for example, 
using Double Frequency Sweep (DFS) [22.23J and Fast Amplitude Modula- 
tion [2^f25] . There have also been improvements made in sensitivity based 
around the inclusion of signal intensity from additional coherence transfer 
pathways [2l?ll2~T] . The Z-filter [28] method ensures that amplitudes for echo 
and anti-echo pathways are co-added with equal intensity under States [25] 
acquisition, providing after phase correction a purely absorptive 2-D spectra. 
Given these improvements, particularly the latter, it is reasonable to assume 
that lineshapes for individual crystallite orientations may be described via 
traditional linear response theory [30, 31~f32"] . Further, allowing for the possi- 
bility of contributions from both homogenous and inhomogeneous broadening 
processes, a complete model includes a linear combination of Lorentzian and 
Gaussian absorption lineshapes with broadening factors A2, Ai: 
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where e < 1 is a free parameter, describing the relative fraction of different 
lineshape functions. It is assumed for the remainder of this work that at- 
tention is restricted to symmetric transitions (eg., 3QMAS experiments) and 
thus devoid of first order quadrupole effects, or that first order effects are ab- 
sent from satellite transitions, the latter ensured by using an accurately set 
magic angle. Finally, it is assumed that experiments are conducted using a 
rotor-synchronized Fl dimension to eliminate spinning sidebands in this di- 
mension [33]. Under these assumptions, the indirect 2irfl m = oj^J and directly 

detected frequencies 27r/2 m = lo^\ have the general fornH] 

<*% = (r " cHC - r, c) (i±^) 

+ A( 4 )(J,r,c)/( % ,a,/3)}, (3) 
where: 

A (0) (I, r, c) = I( I + 1) - 3(r 2 + rc + c 2 ) 

A {A) (I, r, c) = 18/(7 + 1) - 34(r 2 + rc + c 2 ) - 5 (4) 

are spin (J) and quantum transition (r, c) dependent constants. The isotropic, 
second order quadrupole shift t^o is given by the second set of terms in equa- 
tion 3 divided by the Larmor frequency u>o, and contains the quadrupole cou- 
pling constant implicitly: 
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Cq T (5) 



The second order, quadrupolar line broadening is described by function 
f(rj q ,a,/3) [311 • This term is a function of the asymmetry parameter r] q and 
powder angles a, /3, the latter describing the orientation between the Principal 
Axis System (PAS) of the electric gield gradient tensor and the rotor fixed 
frame: 

f(rj q , a, (3) = t(~ 54 ~ 3l lq + 60% cos 2a - 35r/J cos 4a) 

+ (540 + 30r/ 2 - 480r/ g cos 2a + 70r/ 2 cos 4a) cos 2 (3 

(-630 - 35r] 2 + 420r, 9 cos 2a - 35r, 2 cos 4a) cos 4 0\ (6) 



2 Frequency transitions are labeled by r and c. Owing to the dipole selection rule, 
directly detected frequency transitions are always such that r — c = ±1 eg., the 
central transition (—1/2 «-> 1/2). Multiple quantum transitions are such that r — 
c 7^ ±1 eg., the triple quantum transition (—3/2 <-> 3/2). The particular multiple 
quantum transition(s) correlated with the central transition in the course of an 
experiment are determined by the phase cycle. 
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This quantity is a direct consequence of the transformation between principle 
axis frame of the crystallite and rotor fixed frame, in terms of Wigner rotation 
matrices. Assuming experiments are conducted in the fast MAS limit, where 
attention may be restricted to the centerband, a third angle 7 describing the 
rotor orientation with respect to the static field is unnecessary, since the static 
field represents a symmetry axis for spins. Throughout the course of an MQ- 
MAS or STMAS experiment, or via subsequent data processing, the indirect 
dimension frequency fl m becomes fV m = fl m — k x f2 m , the shearing factor k 
often chosen to eliminate the anisotropic frequency component and thus create 
a fully isotropic frequency dimension. The resultant frequency fl' m as well as 
the accompanying bandwidth may be rescaled by a factor 1/(1 + k), according 
to one convention. For ease of comparing spectra arising from different multi- 
ple quantum experiments, this work follows the unsealed representation [35] . 
Note that regardless of the convention followed in presentation and analy- 
sis of spectra, the isotropic chemical shifts ultimately deduced are identical. 
Equation 2 is germane to a single crystallite orientation, with a particular 
isotropic chemical shift, asymmetry parameter and quadrupole coupling con- 
stant. A more general lineshape intensity function for a powdered solid must 
be weighted by crystallite angle distribution G(a,j3). In addition, in the pres- 
ence of disorder, the experimental lineshape is averaged due to distributed 
values of 5 l c s s °,C q ,r] q , described by probability density P(5 l c s s °,C q ,rj q ): 



where M is the total number of chemical sites and A% the individual site 
amplitude. There are two basic aspects to a numerical evaluation of this five 
dimensional integral, including powder averaging over the crystallite orienta- 
tions. In addition, contributions to the overall spectrum from random variates 
C q , 5 l c s s °, and r\ q are weighted by a multi-variate probability distribution func- 
tion, distinct for each site. The former aspect, powder averaging in magnetic 
resonance, is an example of a problem in broader quantum mechanics, eval- 
uating integrals over the unit sphere [5oTI3T] . There exist several reviews in 
the literature with regard to powder averaging in magnetic resonance [SHIES]- 
It is assumed here that the equally probable crystallite orientations within a 
powder have been equally irradiated, and the integral over angles is replaced 
by a sum: 



with various choices for weights Wk and angles a, (3. Under this assumption, 
the contribution of a particular crystallite orientation to the overall inten- 
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sity is proportional to dad(3 sin j3 ie., G(a,/3) = sin/3. The particular powder 
integration scheme used within this work corresponds to the Zaremba-Conroy- 
Wolfsberg (Z-C-W) method [4"Uf41f42| . where angles and weights are chosen 
according to: 

_ 2w(kM a mod Ng) 

" fc ~~ N a 

fa = arccos (l - (9) 
w k = 1 



with N a and M a chosen to satisfy M a = F(n) and N a = F(n + 2), where 
F(n) is the nth Fibonacci number, and index k = 0, 1, N a — 1 . This par- 
ticular approach is considered preferable under fast MAS conditions |39j and 
demonstrates very good convergence versus order n. 

The second major aspect to evaluating equation 7 involves averaging over 
isotropic chemical shift and quadrupole parameters, accomplished via Monte 
Carlo simulation. In general, statistical distributions may be symmetric or 
asymmetric. The nature of the model distribution used in the simulation is 
directly related to the underlying chemical and/or structural disorder. Tradi- 
tional random number generators which create variates according to proba- 
bility distributions are usually one of two types. They may be of the accep- 
tance/rejection type, or rely on transformations of the uniform distribution, 
eg., the Box-Muller method for normal-distributed variables [13]. The latter 
was used here for ease of adaptation to a parallel programming environment. 
By creating Gaussian distributed variates, the integral of eq. 7 over the proba- 
bility distribution may be converted to a summation, and the powder-averaged 
kernel F(fl,f2) simply evaluated as a function of the variates. By the law 
of large numbers, Monte Carlo approximations converge to the true value in 
the limit as the samples N approach infinity. In reality, convergence is slow, 
and the error in using pseudo random numbers is 0(N~ 1 ^ 2 ). This situation is 
improved via using quasi-random numbers such as the Sobol sequence, which 
have an error 0((logiV) m A^~ 1 ) for m dimensions |44j . For the purposes of this 
work, attention is restricted to the bi-variate (m=2) Gaussian distributions in 
5™° and C q (whose random variates are represented by x and y respectively): 

CT g ^°y al 

P(x,y) = e 2 ~^ r ) (10) 



For each chemical site, this distribution is parameterized by site-specific values 
for fi x , a x , fiy, <jy, p, where p is the correlation coefficient between chemical 
shift and quadrupole coupling constant only. At this stage, single values for 
r] q were deemed sufficient to model lineshapes. This was due to an observed 
insensitivity of lineshape simulation to a range of values for rj q . To summarize 
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thus far then, each chemical site % is modeled using ten free parameters a 1 : 



{A}, A|, e\ fj.1 oi, ^, p\ rf Q , A% i = l,,M (11) 



The integral of (7) is replaced by a double summation, in performing powder 
angle and parameter averaging tasks: 



1 



M N N a —1 



N (N _ n EAE E (12) 

where N is the total number of variates x, y for each chemical site i. These 
variates are sampled from a bi-variate Gaussian distribution, using the Box- 
Muller transformation of (Sobol) quasi-random numbers on [0,1). The powder 
angles a, (3 are chosen according to the Z-C-W scheme, as is the number of 
summands (N a — 1). 



2.2 Optimization 



Using the theory outlined thus far, an experimental spectrum may be simu- 
lated and attempts made to optimize the simulation parameters. In reality, 
modeling the underlying parameter distributions implies that at least two 
chemical sites are used in the optimization. Figure 1 is a plot of the RSS func- 
tion obtained by varying only chemical shifts in an optimization for a two site 
MQMAS spectrum. 




Fig. 1. RSS function, sum of squared difference between simulated and experimental 
MQMAS spectrum, as a function of the two isotropic chemical shifts. 

The surface is highly non-convex; the global minima is toward the center of the 
plot, within a larger area containing local minima. Simulated annealing [15] is 
a stochastic method for global optimization highly suited to non-convex RSS 
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functions. The method is analogous to the metallurgical process of annealing. 
The application to the current problem ensures that the iterative procedure 
avoids being trapped within local minima; the overall algorithm applied here 
is as follows: 

(1) RSS function or generalized energy generation, the trace of the Gram- 



E = Trace {(A — B) x (A - B) T } 

where A — B is a matrix of residuals, the difference between simulated 
A and experimental absorption spectra B. If this is the initial step, a 
generalized temperature is defined T « Eq 

(2) Each unconstrained parameter a* is changed by a random amount ±pAa\ 
p sampled from the uniform distribution [0, 1). The corresponding energy 
Ef is calculated as before. 

(3) If Ef < Eq, the change is accepted, else, 

(4) Parameter changes are accepted or rejected in the traditional Metropo- 
lis jl6] scheme, using the probabilistic factor: e~( Ef ~ Eo ^' T 

(5) The process is repeated and the temperature lowered according to some 
schedule, until such time as convergence is reached. 

Implicit to the algorithm is the need to choose an appropriate maximum step 
size Aa* and annealing schedule. To ensure adequate search of the parameter 
space, Aa' was fixed between one and two percent of the starting parameter 
values. The annealing schedule is more subjective and best determined via 
experiment. A common method involves reducing the temperature at every 
step by some amount 5: 



which requires the tuning of 5. Significant gains are made during the early 
stages of the algorithm, during which there is a non-zero probability for energy 
to increase. In order to exploit this feature, 5 was set to approximately 0.5 
and the schedule of equation 13 was re-set every k steps to the current best 
value of energy, a process of rapid annealing and re-annealing. 

2.3 Error Estimation 

In order to give confidence intervals for the free parameters listed in eqn. 
11 optimized in the simulation, strictly speaking the measurement or MQ- 
MAS experiment in conjunction with simulations ought to be repeated and 
statistics created from fitted data. However, owing to the considerable time 
multiple experiments and simulations requires, a more suitable approach to 



mian: 




(13) 
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error analysis is found in statistical re-sampling [47], such as jackknifing or 
bootstrapping [4"8] . In the original jackknife approach, cf)_j is defined as the 
least squared estimate of parameter when the jth data point of n total is 
removed from the set. Pseudo values are created, 

p. = n( j ) -(n- 1)0^. (14) 
with average P and variance matrix Vp\ 

n 

P = ^ J = n- 1 Y J P 3 (15) 

3=1 

1 



nVp 



J- A — 1 



In the present application, this method implies n + 1 non-linear optimizations 
which is still far too time consuming. Fox et al [20] propose a solution in the 
form of an approximate jackknife, which requires instead a single non-linear 
optimization, via a Taylor expansion of the least squares estimate equation 
for <j)j, assuming it is a stationary point for the sum of the residuals. In this 
method, an estimate of the variance matrix Vj is given by: 

n 

Vj = (Z T Z)-^^Jrj(Z T Zr (17) 

3=1 

where: 

Zj = V/(x„0) = | )...^_ /(Xj)0 )|^ _ (is) 

Z T = (z u ...,z n ) (19) 

and Tj is the vector of residuals. The model as presented here consists of ten 
free parameters per chemical site (ie., /=10), so in the case of M chemical 
sites, this corresponds to the creation of a 10M x 10M variance matrix. This 
matrix is evaluated at best-fit parameters a 1 , using the partial derivatives of 
equation 7, listed in appendix A and evaluated as before via summation. 



3 Implementation 



The aforementioned theory was implemented in C, using a number of func- 
tions from the GNU Scientific Library (GSL), as well as the math and standard 
libraries. A single application was written which performs calculations of fre- 
quency equation 3, for each dimension. Further, a multiple k of the direct 
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dimension frequency /2 m is subtracted from the indirect dimension frequency 
fl m , according to the function of the shearing transformation. As stated pre- 
viously, the shear factor k is often chosen to produce a fully isotropic frequency 
dimension and for the examples given here (spin 5/2 and 3QMAS experimen- 
tal conditions) corresponded to a numerical value of 19/12. For each frequency 
dimension, Sobol sequences are generated and used to create bi-variate distri- 
butions of isotropic chemical shift and quadrupole coupling constant according 
to the Box-Muller algorithm. Powder angles are generated according to the 
Z-C-W algorithm. Finally, summation over powder angles, variates and chem- 
ical sites are performed using equation 12. A single OpenMP pragma was used 
to parallelize inner frequency loops, 

#pragma omp parallel for private (h,i) 

using the private declaration on loop indices to prevent a race condition oc- 
curring between separate threads. The OpenMP application programming in- 
terface is essentially a set of libraries and associated compiler directives which 
permits shared memory processing (SMP) on machines with the appropriate 
hardware. The C source was compiled using the GNU C compiler, linking the 
appropriate libraries: 

gcc -04 -o mqmas_opt mqmas_opt . c -lm -lgslcblas -lgsl -fopenmp 

In order to perform optimization of the simulation parameters, the simulated 
annealing algorithm was implemented in an OCTAVE script, mqmasOpt .m. 
This allowed for tuning of heuristic parameters, particularly the annealing 
schedule and size of random fluctuations taken by individual parameters per 
iteration (set to between 2-3% of initial parameter magnitudes). In addition, 
parameter values corresponding to the lowest energy obtained are stored every 
iteration and used for occasional resets. 

It is anticipated that the number of crystallite orientations required for ad- 
equate convergence in a particular simulation will increase with linewidth, 
which in turn is proportional to the quadrupole coupling constant. Fitting to 
a crystalline model compound provides a good means of determining the min- 
imum number of crystallite orientations required for a comparable linewidth. 
Convergence or lack thereof is more easily observed in a crystalline system 
as compared to a more disordered material, which is devoid of the character- 
istic features. In order to test convergence of the powder averaging step for 
a material of interest, a 27 A1 3QMAS spectrum of large-pore aluminophos- 
phate VP 1-5 was acquired using a 11. 7T spectrometer, figure 2. To simulate 
the full 3QMAS spectrum without visible irregularities, 1597 angle pairs (Fn) 
were minimal for quadrupole coupling constants in the range less than 4 MHz, 
as exhibited by the model compound VP 1-5. The Second Order Quadrupole 
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Fig. 2. Spectra are referenced to AICI3 and scaled according to the second conven- 
tion of Amoureux et al [39], where indirect dimension bandwidth reflects spinning 
speed (a) 27 A1 3QMAS VPI-5 spectrum obtained at 11. 7T, using a Chemagnetics 
spectrometer and z-filter/States sequence |28j ; spinning speed 10kHz, bandwidth 
10kHz in each dimension, 64x1024 total points in Fl and F2 respectively, (b) Sim- 
ulation of experimental 3QMAS spectrum (c)Trace along frequency fl = 58 ppm 
showing experimental (lower) and simulated (upper) spectrum (d) Trace along fre- 
quency fl = 53 ppm showing experimental (lower) and simulated (upper) spectrum 
(e) Integrated intensity along isotropic dimension showing experimental (lower) and 
simulated (upper) spectrum. 

Effect (SOQE) parameter^ as determined from the simulation for the tetra- 
3 SOQE = C q \/l + f 
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Table 1 

Results for simulation of hydrous albite 3QMAS spectrum 
Site # 5f so {Hz/ppm) C q (MHz) rj q Rel. Popu lation 

// a [i a 

1 8035/61.7 172/1.3 3.8 1.5 0.21 0.34 

2 8623/66.2 167/1.3 2.9 0.8 0.51 0.66 
Table 2 

Jackknife parameter error estimates for simulation of hydrous albite 3QMAS spec- 
trum 

Site # 5f so (%) C q (%) n q (%) Rel. Popu lation (%) 

1 1.1 4.8 8.1 7.8 8.6 18.0 

2 0.2 1.9 1.4 1.7 0.2 4.0 

hedral region of VPI-5 were 2.8 and 1.3 MHz, which compare favorably with 
literature values [50] . 

Using the same number of crystallite angles, optimized simulations were per- 
formed for the tetrahedral region within a hydrated albite sample, using 200 
quasi-random samples for each of two chemical sites, drawn from two bi-variate 
distributions. Results are displayed in figure 3 and table 1. 

The Gaussian/Lorentzian ratio, correlation coefficient and broadening con- 
stants in both dimensions were constrained to 0.5, 0, and 100 Hz respectively 
and 1000 simulated annealing iterations were performed. The experimental 
spectrum displays regions of both order (narrow, horizontal peaks) and disor- 
der (broad, indistinct). In order to test the validity of the simulated, optimized 
model, jackknife parameter error estimates were determined and are presented 
in table 2. 



The chemical site with narrow distributions (assigned here to crystalline al- 
bite, in good agreement with prior investigations by other workers [51]) has 
corresponding parameters with least error. This may be attributed to a num- 
ber of factors, in this case most likely to the lower signal to noise ratio of 
the disordered region, assigned here to amorphous albite glass. For chemical 
sites with larger quadrupole coupling constants, there is also the possibility 
that due to experimental excitation deficiency, the second order perturbation 
frequency expression breaks down. Finally, the assumptions of a Gaussian sta- 
tistical model may be inappropriate for the system in question. As mentioned 
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Fig. 3. Spectra are referenced and scaled as previously; (a) 27 A1 3QMAS hydrated 
albite spectrum obtained at 11. 7T, using a Chemagnetics spectrometer and z-fil- 
ter/States sequence; spinning speed 10kHz, bandwidth 10kHz in each dimension, 
64x1024 total points in Fl and F2 respectively, (b) Simulation of experimental 
3QMAS spectrum (c)Trace along frequency fl = 96 ppm showing experimental 
(lower) and simulated (upper) spectrum (d) Trace along frequency fl = 91 ppm 
showing experimental (lower) and simulated (upper) spectrum (e) Integrated inten- 
sity along isotropic dimension showing experimental (lower) and simulated (upper) 
spectrum. 

earlier, model distributions reflect the underlying stochastic nature of bonding 
in a disordered material. There is significant evidence [52 53.54J to suggest that 
a more general electric field gradient model for disordered systems as probed 
by 27 A1 NMR is given by the Czjzek model [19]. Future work will be devoted 
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to the incorporation of models such as these into the approach outlined here. 



4 Summary 

Theory has been outlined and an application implemented in the C program- 
ming language that permits the simulation of an MQMAS spectrum, as a 
function of underlying parameter distributions. This simulation relies on the 
use of quasi-Monte Carlo variates to promote convergence and utilizes the 
OpenMP library to permit execution on SMP machines. Owing to the man- 
ner in which random variates are created in the application, the program is 
amenable to High Throughput Computing (HTC) platforms such as Condor 
or PBS. In addition, an OCTAVE script implementing a simulated annealing 
algorithm is used to optimize the simulation, providing reliable estimates of 
NMR parameters. Finally, theory was outlined and implemented for providing 
parameter variance estimates using a jackknife approach. In conjunction with 
the MQMAS experiment, the application described herein enables the char- 
acterization of materials which may vary greatly in the degree of underlying 
chemical and structural order. 
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Appendix A 



Define: 



clb { 



1(1 + 1) + 3/4 



clbi 



18/(7 + 1) + 34/4 + 5 
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db 2 = (r - c)(I(I + 1) - 3(r 2 + rc + c 2 )) 

c /6 3 = ( r _ c ) (18/(7 + 1) - 34(r 2 + rc + c 2 ) - 5) (20) 

then 

dl = dl dfl m | dl df2 m 
di] q dfl m dri q df2 m dr] g 

where: 

df2 m 

dr) q 

clbxy 2 
15120/ / 2 (27- l) 1 

+ cos 4 /3 (-70.0 cos (4 a) rj q - 70.0 rj q + 420 cos (2 a)) - 70.0 cos(4.0a)^ 
-6.0^60.0cos( 2a )}-^^ T? (22) 

dfl m 



{cos 2 /3 (140 cos(4.0a)r? g + 6O.O77, - 480 cos(2a)) 



dr] q 

15120/o/^(2/ - l) 2 { C ° S2 ^ (M0 cos ( 4 - 0a )^ + 60 - ^ " 480 cos ( 2a )) 
+ cos 4 /3 (-70.0 cos (4 a) rj q - 70.0 ^ + 420 cos (2 a)) - 70.0 cos(4.0a)^ 

-6.H + 60.0cos(2a)} + 5/ ^y; )2 



2 



" 15120/!%/ -1)2 { C ° S2 ^ (14 ° C ° S (40a) ^ + 6 °- H " 480 C ° S(2a)) 
+ cos 4 (-70.0 cos (4 a) rj q - 70.0 rj q + 420 cos (2 a)) - 70.0 cos (4.0 a) r, q 

-6.0 Vq + 60.0 cos(2a)} + ^jf^ 1)2 (23) 

where k is the shear factor of the MQMAS or STMAS experiment, 
dl 



df2 rn 

( 

AiPi(x,y) 

V 



(/2 -/2 m )e ^ ^ e 

2ttAiA| 
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2(/2-/2 m )A 1 A 2 (l-e) 



(A? + (/l-/l m ) 2 ) (A 2 + (/2-/2 m ) 2 )' 
dl 



(24) 



Oflr 



( 



AiPi(x,y) 



(/l - fl m )e 



2X 



i e 



2ttA 2 A? 

2(/l-/l m )A 1 A 2 (l-e) 
(A| + (/2 - f2 m f) (A 2 + (/l-/l m ) 2 ) : 



(25) 



sa; 



AiPi(x,y) 



/ (/2-/2 m ) 2 (/l-/l m ) 
m m 



(/2-/2 m ) 2 (/l-/lm) 2 

i e , (/2-/2 m ) 2 e "~^^ e 



V 



2vrAiA| 



Ai(l-e) 



+ 



(A 2 + (/l - fl m f) (A 2 + (/2 - /2 m ) 2 ) 



2AiA^ (1 - e) 



(A? + (/l - fl m f) (A| + (/2 - /2 m ) 2 ) j 



(26) 



^aT 



AiPi(x,y) 



/ (/2-/2 m ) 2 (/l-/lm) 
T—^ T—^ 



e 2A 2 



V 



2vrA 2 A 2 
A 2 (1 - e) 



_(/2-/2 m ) 2 _(/l-/l m ) 2 

^6 | (/2-/2 m ) 2 e ~^~e 



27rA 2 Af 



(A 2 + (/l - /l m ) 2 ) (A 2 + (/2 - /2 m ) 2 ) 



2A 2 A 2 (1 - e) 



(A 2 + (/l-/l m ) 2 ) (a 2 + (/2 - /2 m ) 2 )' 



(27) 



di 
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A i P i (x,y)F i (fl,f2)x 

( p p{x-p x f (x-n x )(y-n v ) 2p 2 (x-p x ) p(y - p y f \ 



V(l-p 2 ) (l-p 2 )V x (l-p*)a x a y (l-p*)a x a y (l-p 2 )%, 

dl 

da y 

AiF t (fl, f2)P i (x, y) (2p{x- p x ) (y - p y ) 2 (y - p y f 



2(1-P 2 ) V °l 
A i F i (flJ2)P i (x,y) 



a y 

dl 

dp y 

_ A l F, l (flJ2)P t (x,y) ( 2p (x-p x ) _ 2 (y - p y ) \ 

2(1"P 2 ) V °x<Ty °l J 

dl 

da x 

A^fl, f2)Pi(x, y) (2p(y- p y ) (x - p x ) 2(x- p x ) 



2(1-P 2 ) V 
A i F i (fl,f2)P i (x,y) 

dl 

A i F i {flJ2)P i {x,y) (2p{y-p y ) 2 {x - p x )\ 



2(l-p 2 ) \ a y a x 



(29) 



(30) 



(31) 



(32) 



dl 

— = F i (flJ2)P i (x,y) (33) 
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